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The abiUty to generate complete, or almost complete, chaotic mixing is of great interest in nu- 
merous applications, particularly for microfluidics. For this purpose, we propose a strategy that 
allows us to quickly target the parameter values at which complete mixing occurs. The technique 
is applied to a time periodic, two-dimensional electro-osmotic flow with spatially and temporally 
varying Helmoltz-Smoluchowski slip boundary conditions. The strategy consists of following the 
linear stability of some key periodic pathlines in parameter space (i.e., amplitude and frequency of 
the forcing), particularly through the bifurcation points at which such pathlines become unstable. 

PACS numbers: 47.51. -fa, 47.52. 47.61. Ne 



I. INTRODUCTION 



As it is well-known, the generation of efficient and complete mixing is particularly challenging for Stokes flows, 
for which the Reynolds number is very low {Re < 1) and turbulence does not exist. These include large scale flows 
for fluids of high and/or low velocity as well as small scale flows [T]. At miniature scale, efficient mixing leads to 
important applications, particularly lab-on-chip devices that can incorporate one or several laboratory functions on 
a small surface (with a size of approximately one millimeter square or less). In practice, such devices are often used 
as biochemical reactors for which complete mixing is essential in order to bring different reagents in contact of one 
another and trigger chemical reactions between them. 

For the past three decades, the kinematics viewpoint of fluid mixing, that is the behavior of trajectories of passive 
or advected particles, has received much attention (see ^2J). This is through this viewpoint that the intimate relation 
between fluid mixing and chaotic advection, also referred to as chaotic mixing, has been explored in detail. Specifically, 
it is now well-understood that chaotic mixing occurs through successive mechanisms of stretching and folding of 
material lines and that it increases exponentially the contact area between the different fluids or reagents to be mixed. 

The generation of chaotic advection in a small Reynolds number flow, particularly in a microfluidic device, is 
generally achieved by adding a degree of freedom to a two-dimensional incompressible base flow. Such a degree of 
freedom takes the form of time dependence [3HB] or a third dimension [71 |H] . 

Generally, this time dependence or third dimension is periodically introduced into the system via passive techniques 
based on altering the device geometry [SHU], active techniques based on forcing (e.g., [T^HTT] ). or the combination of 
both [T5H^ . In this work, we use the same base flow as in Ref. [13] with a time periodic dependence introduced into 
the system via forcing which is created by temporally varying the slip boundary conditions along the channel walls. 
While the forcing was created by using a discontinuous function in the form of a switch in Ref. 13J, our forcing is 
created by a smooth time-periodic oscillation of the slip boundary conditions. The smoothness of the system is easier 
to handle computationally and experimentally. 

In this paper, we focus on creating and spreading chaotic mixing within the entire channel. Our goal heads towards 
the opposite direction of other works that target the total elimination of chaotic mixing or those that target partial 
mixing at a desired location and size in phase space [H]. Our strategy follows the one in Ref. [23] j where the global 
extension of chaotic mixing can be summed up by the linear stability of a few key periodic pathlines. With this 
strategy complete or at least almost complete chaotic mixing is achieved. 

The paper is organized as follows. The physical model as well as its assumptions and corresponding dynamical 
system are described in Sec. [llj Section III shows the complex relationship between the system parameters (i.e., 
amplitude and frequency of the perturbation) and the mixing characteristics. In particular, this section outlines 
how the trial and error methods commonly used to flnd complete chaotic mixing is burdensome and computationally 
demanding. Indeed, for only one chosen pair of parameter values, such an approach (based on Poincare sections, 
flnite time Lyapunov exponent maps and mixing indices) requires the integration of the dynamical system for a large 
number of initial conditions over long periods of time. In Sec. IV we present a strategy capable of determining the 
set of parameter values leading to complete, or almost complete, chaotic mixing in a much more efScient and less 
time consuming fashion. In contrast to the previous methods our method focuses on a few invariant structures, i.e., 
key periodic pathlines over a very short period of time. It allows us to identify the parameter values for which well 
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FIG. 1: Schematic of the microfluidic channel subjected to a steady uniform electric field Eq in the downstream direction. In 
addition, the channel is equipped with a series of electrodes embedded within its walls. The channel can thus be decomposed 
into a series of generic periodic cells defined by a set of four adjacent electrodes. The black dashed line represents the boundary 
of a generic cell. 

spread chaotic mixing takes place. 

II. PHYSICAL MODEL 
A. Assumptions 

We consider a straight two-dimensional microfluidic channel of height 2H, filled with a weak electrolyte solution 
of electrical permittivity e^i and kinematic viscosity v. Electroosmotic pumping of the fluid is achieved by applying 
a steady, uniform electric field Eq in the streamwise direction by means of two electrodes placed upstream and 
downstream of the channel (see Fig. [T]). When the electrodes are energized, it is known that an electroosmotic flow 
is generated due to the accumulation of counterions near the charged channel walls and their downstream migration. 
While the electrical double layer containing the counterions is typically thin compared to the channel height, it 
nevertheless drags the flow with it, thus resulting in a quasi-uniform flow velocity profile. Assuming that the electrical 
double layer is indeed sufficiently small, the magnitude of the flow velocity is defined by the slip velocity Vhs at the 
walls and its expression is given by the Helmoltz-Smoluchowski formula 

Vhs = ^, (1) 

V 

where C is the potential difference across the electrical double layer, also called the ^ potential. It is clear from 
the above expression that the flow velocity, in both amplitude and direction, can be controlled by adjusting the C 
potential. As in Ref. this is performed by placing a series of adjacent electrodes of length 2H within the upper 
and lower walls f2T (see Fig. [ij. By varying the voltage applied to the side electrodes, the slip velocity Vhs at the 
walls can be varied in time, independently of -Eg [25n28i . 

Due to the small scale of the microfluidic device (with a height typically of the order of 10 — 100 micrometers) 
and the relatively weak electroosmotic velocity (of the order of 10~^ m.s~^), the Reynolds number is typically much 
smaller than one. It is thus reasonable to model the flow inside the channel as an incompressible Stokes flow. The 
boundary conditions, which consist of a periodic slip velocity (in both space and time) at the upper and lower walls 
{y — ±H, respectively), are expressed as 

Vx (x, ±H, t) = Vhs {±Ra2H (x) + e cos ojt) , (2) 

where Ra2H [x] — A/n sin ((2'T' — 1) x/2H) / {2n — 1) denotes the periodic square wave function of period 2H 

bounded by the values (—1, 1) and A = {e,a;} stands for the amplitude and frequency of the forcing. The unsteady 
part of the slip boundary conditions is thus characterized by two parameters, the amplitude e and the frequency lo. 
As we show below, control of chaotic mixing will be achieved by adjusting these parameters. 
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B. Velocity field: Dynamical system 



Due to the linearity of the Stokes flow problem, the velocity field Vx can be decomposed into a steady component 
Vo obtained for e = (see Fig.[2]-a), and an unsteady one, Vi, so that 

Vx {x, y, t) = Vo {x, y) + Vi {x, y, t) , (3) 

where the variables have been made dimensionless by using Vhs ^nd H/Vhs ^ length, velocity and time scales. 
In the steady case, one obtains recirculating rolls in between the upper and lower walls (see Fig. [2]-a), as shown by 
[24| . The unsteady, perturbed case is then generated by means of the unsteady slip velocity obtained by temporally 
varying the C potential at the walls. 

For an incompressible flow, the dynamics of a passive particle can be studied by using the stream function formu- 
lation. Its velocity is expressed as 

^^.^_5V^x^ (4) 
dy 

dipx{x,y,t) 

""^y^^x — ' 

with 

'4^\{x,y,t)^^o{x,y) + i)i{x,y,t). (6) 

This system is Hamiltonian, with the stream function t/'a acting as the time dependent Hamiltonian and the physical 
coordinates {x, y) playing the role of conjugate variables. The stream function is composed of a steady component, 
■00 (corresponding to Vq), and an unsteady one, ■01 (corresponding to Vi). The steady part is determined by solving 
the steady Stokes equation with periodic slip velocity at the walls. Such a problem can also be approached by means 
of the stream function formulation, i.e., by solving the biharmonic equation 

VVo = 0, (7) 

where the stream function of the steady case, tp^, satisfies impermeable and periodic slip conditions at the walls 

^Pq{x,±1) = 0, (8) 

^(x,±l) = ±Ra2H{.x). (9) 

The biharmonic equation ([t]) is then solved by the method of separation of variables, the solution taking the form 

oo 

00 {x, y) = ^ a„X„ {x) Yn (y) , (10) 



where 



X„(x) = sinf^), (11) 



V 2 



tanh 



((?!i_i)£) cosh (fili^) - ysinh 



(2n— l)7ry 
4 



^"(y) = ^ ' ^ (12) 



cosh ' (2^^ 



cos (riTr)) cosh^ 



;^)((^)+cosh2(^)tanh(^)) 



(13) 



For further details on the derivation of the solution ipo, we refer the reader to reference [24] . 

The unsteady component (or perturbation), generated by the oscillating slip conditions of amplitude and frequency 
A — {e, uj} at the walls, act in the y direction only and can be trivially derived as 

%pi{x,y,t) ^ eycosujt. (14) 



The resulting velocity field in a generic cell is represented in Fig. [3^-c in space and time. 
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FIG. 2: a) Velocity field and contour levels of the velocity intensity inside a generic cell for the steady case (e — 0). b) 
Streamlines of the flow, highlighting the heteroclinic pathlines (dashed lines), the hyperbolic fixed points (crosses) and the 
elliptic fixed points (thick dots). 



III. MIXING CHARACTERISTICS AS A FUNCTION OF PARAMETERS 



A. Integrable case 



The steady case, corresponding to e = when the time dependent voltage applied to the wall electrodes is turned 
off (or equivalently, without time periodic slip conditions), is characterized by one invariant of the dynamics, the 
stream function V'o G [-i>7nax,i>max] given by Eqs. ([To|-([l3l). 

In this case, the flow consists of a series of periodic cells, each cell consisting of two side-by-side recirculating rolls 
rotating in opposite directions as represented in Fig. ^p. Within a generic cell, pathlines are organized into two sets 
of closed lines around elliptic fixed points, r±^^^^ located at (1,0) and (3,0). In addition, heteroclinic pathlines, 
Tq, connect hyperbolic fixed points located at the corners of a generic half cell i.e., (0, —1), (0, 1), (2, —1), (2, 1) (see 
Fig.[2}b). 

For 'i/'o 7^ the corresponding pathlines are closed lines of constant "00 7 denoted by F^q. Following the non-mixing 
structure of the pathlines, a small amount of dye injected into the flow evolves into a circular-like pattern without 
spreading throughout the channel. Clearly, no mixing is generated in this case. 



B. Slightly perturbed case 



With a time dependent perturbation introduced by the oscillatory slip conditions, the system acquires an additional 
degree of freedom, thus making the generation of chaotic mixing possible. We first consider the weakly perturbed 
system, that is the flow subjected to a perturbation of small amplitude (e <^ 1). Figurejsji displays a Poincare section 
of the dynamics given by Eqs. Q-llsl), that is a stroboscopic map of period 27r/a; modulo the length of the cell. The 




FIG. 3: a-c) Velocity field and contour levels of the velocity intensity inside a generic cell for the unsteady case with e = 0.10 
and uj = 0.90 at time t = 7r/2aj, nuj, Stv/Auj. d) Corresponding Poincare section. 
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phase space in each generic cell consists of a small chaotic mixing region created by the destruction of the heteroclinic 
pathlines present in the unperturbed system. Deeper within the center of the two half cells, two symmetric sets of 
quasi-periodic pathlines revolve around elliptic periodic pathlines (labeled as O^^^^ and 0_^^^^). In such regular 
islands no mixing occurs. 



C. Qualitative observation of the chaotic mixing 

The perturbation -01 generated by the time periodic slip conditions at the walls plays a crucial role as it is the part 
responsible for the creation of chaotic mixing. Indeed, such a periodic perturbation adds a new degree of freedom to the 
system under the form of time t, thus making chaos possible. Chaotic zones in the phase space, which here corresponds 
to the real space, are created by a web of higher order resonances between the base flow ipo and the perturbation 
ipi. The structure of the phase space is then composed of both chaotic and non chaotic zones intertwined in a fractal 
structure, with a multitude of non- mixing islands as various scales, as shown by Poincare sections (see Fig.[4^-c). The 
presence of large non-mixing islands prevents the spreading of chaotic mixing throughout the entire fluid domain. 

D. Quantification of the intensity and spreading of chaotic mixing 

1. Local quantification: Finite time Lyapunov exponent map 

In order to gain further insight into both the intensity and the spreading of chaotic mixing within the channel, we 
compute the finite-time Lyapunov exponents (FTLE). The technique consists of associating a FTLE C with an initial 
condition Xq — {xq, yo,t = 0). 

First, we consider the time evolution of the Jacobian J* {x,y) given by the tangent fiow and the matrix of variations 
Wx as 

dp 

-j^ = VVx{x,y,t)j\ (15) 
where J'^ = / is the two-dimensional identity matrix. The FTLE map for a finite-time i = r is then defined as 

'C(Xo,r) = -ln|7„a2; (^o)l , (16) 

T 

where 7max is the largest eigenvalue (in norm) of the Jacobian J"^. 

The FTLE map — )■ C (-X'g, t) for some given time t allows us to distinguish between the initial conditions leading 
to the presence of regular (non-mixing) islands, i.e. regions associated with small FTLEs, and the initial conditions 
leading to chaotic mixing characterized by larger FTLEs. The structures in phase space are then easily identified and 
the relative sizes of the regular (non-mixing) islands determined. This tool can be used to not only determine the 
phase space structures but also quantify the degree of the mixing produced. Indeed, large FTLEs indicate a strong 
mechanism of stretching and folding, which is the archetypical mechanism of chaotic mixing. Due to the symmetry 
of the system, the domain of the FTLE map reported below has been chosen as a periodic cell of the channel, i.e. 
{0 < a; < 4, -1 < y < 1}. 

Figure [5] shows the FTLE map at time r = 2007r/w for the set of parameter values e — 0.225,0.975,1.475 for 
Lo = 1.800. It is clear that in Figs. [5^ and [5]: (see also Figs. |4^ and[4j;) chaotic mixing is far from being complete. 
Indeed, around regular islands mixing is weak (i.e., cold colorj, while away from these islands mixing is very strong 
(i.e., hot color). In Fig.[5]D (see also Fig. |4]d), chaotic mixing is clearly well spread, as indicated by the more uniform 
FTLE map across the cell. 



2. Global quantification: box counting method 

Another way to quantify the degree of mixing as a function of spatial location is by determining the mixing index 
M through the box counting method (see Ref. [29]). It offers the advantage of being rather easy to implement, fast 
and rather cheap in computing power. For this, we follow Np advected particles and divide the domain into x Ny 
boxes or cells. At each time, the number of particles Ui is computed in each box, and therefore the fraction of the 
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FIG. 4: Poincare sections obtained for the amplitude e — 0.975 and the frequencies uj — 1.500,1.800,2.250 (a-c) of the 
perturbation for 100 initial conditions taken on y — over an integration time T = 500 x 2tv/(jj. 

total number of particles, or particle rate r^. Given the number of particles rij inside each box i, the computation is 
performed as follows. 

Ti = — if rii < Up, (17) 

Up 

r, = I if rii > Up, 

where Up is the average number of advected particles i.e., Up = Np/ (NxNy). After computing the fraction of particles 
in each box and at each time t, the time evolution of the mixing index M (t) is calculated by taking the average over 
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FIG. 5: Finite time Lyapunov exponent maps for e = 0.225,0.975,1.475 (a-c), and to — 1.800. The integration time is 
r — 2007r/cj over a square grid of 400 x 200 initial conditions. 



all the boxes, i.e., 



^ z— 1 



A mixing index converging towards zero {Yimt^^oo M (t) = 0) indicates an extremely weak mixing process, while a 
mixing index converging towards one (limt_>+oo M (t) = 1) corresponds to a perfect mixing process. 

Figure [6] illustrates the mixing index M versus the time t expressed in terms of the number of periods 27r/cj, for three 
sets of parameter values (e = 0.975; w = 1.800) (diamond), (e = 0.975; w = 2.250) (square) and (e = 0.975; w = 1.800) 
(triangle). In the first case, no non-mixing islands of important size are present (see Fig. [4|d), while in the other 
two cases one observes the existence of important non-mixing islands. From Fig. [6] it is clear that the most efficient 
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FIG. 6: Mixing index M versus time t in terms of the number of periods 2tt/uu, computed by tracing Np = 200, 000 advected 
particles over a cell domain divided into N^c x Ny = 3200 boxes, for the parameter values (e = 0.975; a; = 1.800) (diamond), 
(e = 0.975; cj = 2.250) (square) and (e = 0.975; = 1.500) (triangle). 

mixing (in the sense of a mixing index close to one) is obtained in the first case. In the former case where the chaotic 
mixing is well spread, Fig. [Gj^diamond) shows that the mixing index M quickly reaches a high asymptotic value (after 
about 30 periods), namely M !v 0.95. In Fig. [g] (square), the asymptote is also reached quickly, but its value is much 
lower, indicating only partial chaotic mixing. In Fig. |6] (t riangle), M remains at a very low value, indicating poor 
mixing throughout the channel. 

The previous trial and error methods relying on the Poincare sections, FTLE maps and the mixing index to find 
the parameter values leading to complete mixing are extremely precise and rich in information. Unfortunately, such 
analyses require the integration of numerous initial conditions over long time periods for each pair of parameter values, 
which is extremely computationally demanding. In what follows we propose a much more efficient and refined strategy 
to determine the set of parameter values leading to well spread chaotic mixing. Specifically, instead of investigating 
numerous trajectories over long time periods at arbitrarily chosen parameter values, we focus only on a few specific 
invariant structures of the system, e.g., key periodic pathlines, and track the evolution of their linear stability as a 
function of the parameter values. 

IV. CONTROLLING MIXING WITH KEY PERIODIC PATHLINES 
A. Tracking the stability of key periodic pathlines 

While a chaotic zone usually consists of a chaotic sea containing regular islands surrounding remaining elliptic 
pathlines, we seek to eliminate such non-mixing islands in order to obtain almost complete chaotic mixing. 

Our approach is based on the fact that a variation of the parameters, even small, is capable of influencing the 
dynamics in a significant manner. Our goal is then to identify the range of parameters for which complete, or almost 
complete, chaotic mixing occurs. For this purpose, we study the linear stability of some key periodic pathlines in 
parameter space, with such pathlines being chosen from the slightly perturbed system (e ^ 1). 

Specifically, the key periodic pathlines are selected as the twin pair of periodic pathlines 0±^^^^ located in the 
middle of the two regular islands (see Fig. ^) . These periodic pathlines come from the preservation of the strong 
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FIG. 7: Color on line. Time evolution of a periodic pathline (red) with its projection on the {x, y) plane (black) and deformation 
of a small disc of neighboring initial conditions (blue): a) Stable elliptic periodic pathline obtained for the parameter values 
(e = 0.825, oj = 1.800), showing the confined deformation of the neighboring disc; b) Unstable hyperbolic periodic pathline 
obtained for the parameter values (e — 0.975, tj — 1.800), showing the deformation of the initial neighboring disc consisting of 
successive stretching and folding mechanisms. 



fixed elliptic points under slight perturbation while the pathlines near the fragile heteroclinic pathlines Fg 

have been broken down to become chaotic. 

The stability analysis of a periodic pathline enables us to predict how a small disc of initial conditions around the 
periodic pathline evolves (see Fig. [?]). For a stable periodic pathline the disc is slightly deformed (see Fig. [7^), while 
for an unstable hyperbolic periodic pathline, the disc stretches and compresses exponentially in different directions 
(see Fig. [7]d). As it is well-known, such an unstable periodic pathline can generate chaos by an infinite number of 
interactions between its stable and unstable manifolds. 

The periodic pathline under investigation O — {{x (t) , y (i))}jgjg -^j of period T is numerically tracked in parameter 
space using a Newton- Raphson scheme (see Ref. [30] )■ In order to analyze its linear stability properties we consider the 



time evolution of the Jacobian J* {x,y) given by the tangent flow (see Eq. (15) ). The linear stability of the periodic 
pathline O is then given by the spectrum of the two-dimensional monodromy matrix J-^, where T refers to the period 
of the pathline. Indeed, J^j i^o) — ^gx *^^ \x=Xo describes the deformation at time i = T of an infinitesimal sphere 
of neighboring initial conditions surrounding the periodic pathline O starting at Xq — {xQ,yo,t = 0). Since the fiow 
is incompressible, i.e. det J"^ = 1, these properties can be represented by the trace of the monodromy matrix or, 
equivalently, by the value of the Greene's residue [3T1 [32] given by 

2 - tr 7^ 



In particular, if Rq (A) g]0, 1[ the periodic pathline is elliptic and (in general) surrounded by an elliptic (non-mixing) 
island; if Rq > 1 or Rq < O, the pathline is hyperbolic and thus likely to lead to the appearance of chaotic behavior 
around it; and if Rq = and Rq = 1, it is parabolic. 
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FIG. 8: Residue of the key periodic pathline Oi/,„„^ as a function of the perturbation frequency uj and amplitude e (a-b). 



B. Destabilization of key periodic pathlines 

In what follows we select the two periodic pathlines located at the center of the non-mixing regions in Fig. [2]d and 
[sji as key periodic pathlines. They correspond to the values of the stream function tpQ (in the absence of forcing) equal 
to zbV'max- Hereafter, we refer to such pathlines as 0±^^^^. Since Oi/)„„^ and O-^^^^ are symmetric, we restrict our 
analysis to one of them only, specifically 0^„„^ . 

Figure [8] displays the linear stability property of Ovmoi ^ function of the amplitude e and the frequency lo of the 
perturbation. Figure'^ 
value 1.800. Notice t: 



exhibits the residue value Ro^^^^ for e G ]G, 2[ while the frequency w is held constant at the 
lat the linear stability property depends on the value of e in the following fashion. 



• For e e ]0, ~ 0.85[, the residue value Ro^^^^ G ]0, 1[, indicating that 0^„„^ is elliptic. Consequently, the 
periodic pathline 0,p is surrounded by quasi-periodic pathlines forming a regular (non-mixing) island (see 
Fig.il). 

• For e S ]cehi^he ~ 1-15[, Ro^^^^ < 0; indicating that 0^„„^ is hyperbolic. The destabilization of this 
key periodic pathline creates complete, or almost complete, chaotic mixing (see Figs, [sj^, |4|d). Here the word 
"complete" (or "almost complete" ) is used in the sense that although very small non-mixing islands could appear 
they would be swept away by molecular diffusion due to their small size. 

• For e e ]e^e,2[, 0^i,,„„^ is again elliptic, characterized by surrounding regular (non-mixing) islands (see Fig.jsJ;). 

Figure [sJd shows the linear stability evolution of Oi/,,„„^ for cj e ]1,3.2[ while the amplitude e is held constant at 
the value 0.975. As in the previous case, the linear stability property of O^^^^ follows a similar scenario, as we now 
describe. 



First ^ is elliptic for lu G ]l,io^ii « 1.61[, indicating the presence of two large regular (non-mixing) islands 
(see Fig."^). 

Then, O^^^^ bifurcates from elliptic to hyperbolic at w = iOeh, and bifurcates back from hyperbolic to elliptic at 
Lo = oj/ie ~ 1.93. For w S Jweh, w/,,e[, O^^^^ stays hyperbolic, indicating complete, or at least almost complete, 
chaotic mixing (see Fig. ^ 
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FIG. 9: Color on line. Intersections of the stable manifold, W", (blue) and the unstable manifold, W^, (red) of the pair of twin 
periodic hyperbolic pathlines 0^„„^ for the perturbation amplitude e — 0.975 and frequency to = 1.800. 

• Beyond the bifurcation point ujhe, O^^^^ is again elliptic, indicating that complete chaotic mixing does not take 
place due to the presence of regular (non-mixing) islands (see Fig. |4j;) . 

C. Existence of well spread chaotic mixing 

The existence and spreading of chaotic mixing throughout the channel can be further studied by computing the 
stable and unstable manifolds of the key periodic pathlines O^^^^ in a (unstable) hyberbolic state. Figure [9] shows such 
stable and unstable manifolds and VF" together with their transversal intersections for 0^„„^ at the parameter 
values (e = 0.975, a; = 1.800). The transverse intersections between and imply the existence of a horseshoe 
map (Smale-BirkhofF theorem [33]), thus indicating that the system displays chaotic behavior. The horseshoe map 
can, in turn, be viewed as an archetypical chaotic map [31]. Note that the lobs formed by the intersections between 

and are of significant sizes and are thus likely to lead to a well spread chaotic mixing region. 
The well spread chaotic mixing region is illustrated in Fig. |10[ where the time evolution of two sets of adverted particles 
has been computed. Figure [TO] clearly shows that in less than twenty periods two sets of advected particles, each 
initially concentrated at the center of the two half-cells, eventually become widely mixed and well spread throughout 
the entire cell. 



D. Determination of the chaotic mixing domain in parameter space 



By computing the value of the residue of the periodic pathline O0„„^ as a function of both the amplitude e and 
the frequency oj, we now determine the range of parameters for which the key periodic pathline 0^„„^ is hyperbolic 
and thus unstable. In Fig. 11 the dark area corresponds to the parameter values for which 0^„„^ is hyperbolic, with 
its boundary corresponding to the bifurcation curve through which 0^„„^ switches from being elliptic to becoming 
hyperbolic and vice versa (solid black line) . From this linear stability property of Oi/.,„„^ , we deduce the region of the 
domain which exhibits complete (or almost complete) chaotic mixing. 

We clearly observe that this parameter region of complete (or nearly complete) chaotic mixing stretches in a nearly 
linear fashion (see dashed line in Fig. 111. It thus follows that we can determine a linear relationship between the 



amplitude of the perturbation and its frequency in order to obtain complete chaotic mixing, which we found to be 
e = 0.83a; — 0.53. Since this parameter region stretches rather widely around such a linear median (e.g., up to 37% 
for Ld = 1.6), the previous relation could provide a convenient and robust guide to experimentalists desiring to select 
right away the optimal values of the system parameters. 
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V. CONCLUSION 



In this paper, we have presented an efficient strategy for locating the complete, or almost complete, chaotic regime in 
parameter space, i.e., as the amplitude and frequency of the perturbation is varied. In contrast with the common trial 
and error methods requiring the integration of numerous trajectories over long periods of time, the strategy presented 
here focused on a few invariant trajectories, i.e., key periodic pathlines over a short period of time. Specifically, 
we have studied the evolution in time of the linear stability of these key periodic pathlines as a function of the 
system parameters and predicted the spreading of chaotic mixing as such pathlines destabilize. We have also shown 
both qualitatively (using Poincare sections and stable and unstable manifolds intersections) and quantitatively (using 
a finite time Lyapunov map and a mixing index) that such a local approach leads to an accurate prediction of 
complete (or almost complete) chaotic mixing, while being more efficient than the trial and error methods. Finally, 
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FIG. 11: Bifurcation curve for the key periodic pathline 0^„,^^ (solid black line), displaying the parameter values for which 
O^rna^ Is hyperbollc (and thus unstable), and for which complete chaotic mixing is expected (grey domain). A linear fit of the 
median is also shown (dashed black line). 



for the electroosmotic flow considered in this paper, we have determined the sub-domain of parameter values producing 
complete (or almost complete) chaotic mixing. Given the quasi-linear shape of the sub-domain of parameters for which 
complete (or almost complete) mixing takes place, we have also proposed a linear relation between the parameters 
which should be useful as a guide to experimentalists who can then readily adapt the system parameters to optimal 
values. 
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